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Abstract 


A transient  electromagnetic  field  in  free  space  is 
completely  specified  when  the  initial  values  of  the  electric  and 
magnetic  fields  are  given.  Green's  function  for  the  scalar  wave 
equation  can  then  be  used  to  find  the  field  at  later  times.  A 
group  of  computer  programs  that  implement  these  equations  and 
process  the  output  are  presented  in  this  report. 

Key  words:  computer  programs;  conservation  of  energy  and 

momentum;  el ectromagneti c field  propagation; 
finite  light  pulses;  graphic  displays  of  energy 
density;  transient  electromagnetic  fields. 
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In troducti on 


1 . 

In  order  to  compute  transient  electromagnetic  fields 
scattered  by  conductors  or  dielectrics,  as  required  for  the  Wave 
Optics  program  at  the  National  Bureau  of  Standards,  we  have  to 
determine  the  incident  fields  specified  by  their  initial  values. 
These  incident  fields  are  free  fields  that  obey  Maxwell's 
equation  without  sources  in  all  of  space,  disregarding  the 
scatt erer . 

The  programs  described  in  this  report  allow  us  to  compute 
the  free  fields  at  arbitrary  points  and  times  from  given  initial 
values  of  the  fields. 

A simple  example  is  the  plane  wave:  the  fields  are  constant 

A 

over  planes  perpendi cul ar  to  a direction  of  propagation  n.  The 

4 

initial  electric  field,  is  a vector  field  perpendi cul ar  to 

n that  is  a function  of  a single  variable?  the  free  electric 
field  at  other  times  is  then  given  by 

E(x,t>  = Eq(£  - ct)  (1) 

where  c is  the  speed  of  light  in  vacuum  and 

4 A A -4 

£ = x * n,  n Eq  =0.  (2) 

The  magnetic  field  of  this  plane  wave  is  then  given  by 

cB(x,t)  =nxE(x,t).  (3) 

For  fields  other  than  plane  waves,  the  incident  field  can  be 

4 4 

computed  by  integrations  from  the  initial  values  of  E and  B.  The 
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solution  of  the  free-field  Maxwell  equations  is  shown  in  section 


2,  computer  programs  developed  to  do  the  numerical  integrations 
are  described  in  section  3,  and  a few  sample  outputs  are  shown  in 

section  4.  The  listings  of  the  computer  programs  are  given  in  an 

appendix . 

2.  tLEQe^gat  i^on  gf  Localized  Pulses 

We  have  shown  C 1 3 how  the  retarded  Green  function  for  the 
scalar  wave  equation, 

s‘0> <*,*■)  = *!t  - I 'V  ^'l/C).  (4) 

can  be  used  to  express  the  incident  electromagnetic  fields  in 

terms  of  their  initial  values  given  at  time  t = 0 for  all  of 

space.  The  symbol  x stands  for  the  four-vector  with  components 

•4 

x , m = 0, 1,2,3,  that  is,  x = (ct,x),  and  5(t)  is  the  Dirac  delta 

M 

function.  The  resulting  integrals  are 

-*  f r 2 -*  (0) 

e ( x > =-  | ydV' I (i/c  )eq<x ' )as^  <x,x' ) /at' 

- Bq(x')  X ?'G‘lJ'(x,x')]t,=0  , (5) 


Numbers  in  square  brackets  indicate  the  literature  references  at 
the  end  of  this  report. 


B(x>  = 


- (1/c  >|vdV 


fvdV'[l0(^') 


~ , _iy;  . , . 

v Br  ( x , x ) 


-+  -+  (O)  1 

+ bq(x/ )asR  i*,x' )/at']t,=0  , <6) 

where  v'  is  the  gradient  operator  with  respect  to  the  variable 
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x'.  We  substitute  the  Green  -function  (4)  into  eqs.  (5)  and  (6) 
to  obtain 


-»  r ri  - -»  sf 

E(x)  = JvdV'  [^E0(x'>  — 


(t  - R/c) 
4nR 


+ 


“f  -4 

B ^ ( x ' ) x 
u 


(S' ( t - R/c)  . 5 (t  - R/c ) \ £1 
( 4nRC  4nR2  )RJ 


(7) 


B (x  ) 


(S' (t  - R/c) 
\ 4nRc 


+ 


S (t 


- R/c) 
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4ttR 


R 


- B. 


t - R/c) 

ix  >■ 


4ttR 


(8) 
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where  R = x — x' , R = IRI,  and  R = R/R.  We  integrate  over  R and 
are  left  with  integrals  over  a sphere  of  radius  R = ct  centered 

”4 

at  the  field  point  x,  known  as  the  information-collecting  sphere. 

-4 

We  can  use  the  solid  angle  from  x and  express  the  fields  as 


ct  5 


(9) 


E(x)  = j- 

4rr 


M 


K - + 


“4  ~ T “+  "4  1 

2cB0  ” R-V,cBq)  x r|r 


-r  If  f r -4  -+  -+  -4  -+  -+T 

cB(x)  = 4 — dO  cB  - R-v'cB,  - (2Ert  - R-y'Ert)  x R _ . . (10) 

^nj  LO  0 O O JR  = ct 


The  fields  and  B^  are  arbitrary  solenoidal  fields,  that 
is,  they  satisfy  the  constraint  equations 


EQ(x)  = 0, 


(11) 


v ■ BQ(x)  = 0.  (12) 

A simple  way  to  find  a solenoidal  field  is  to  take  the  curl  of 
another  vector  field. 

We  would  like  the  field  to  be  spatially  localized,  but 
otherwise  to  look  somewhat  like  a plane  wave  propagating  in  the 
2-direction,  We  restrict  the  initial  fields  to  the  region  z < 0 
and  we  assume  that  they  decay  essentially  like  a gaussian  in  the 
x — and  y-di recti ons.  We  choose 

E„(x)  = A expj^  — k (x^  + yJ-)j(yi  — xj)<t>'(z),  (13) 
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cBq<x)  = A exp 


[- 


k (x-^  + y2)j^(xi  + yj)$y(z> 


- 2k  [l  - k ix2  + y2)  J 4>  <z  > y. 


(14) 


where  $ is  a function  that  vanishes  for  z £ 0 and  has  a 


continuous  first  derivative. 


These  expressions  can  be  used  to 


compute  the  dyadics  that  are  the  gradients  of  and  B^. 

The  initial  values  of  the  Roynting  vector  and  energy  density 
are  proportional  to 

Eq  x cBq  = A^exp^  - 2k  (x*"  + - K (x“  + y^)  j <x  i + y j ) 


2 2 2~\ 

+ (x  + y )$>'  k>. 


(15) 


->'?  r 2 2 

E^T  + c^B^  = A^exp  — 2k  (x  + y ) 


K2(* 


2 2 

+ y )<*' 


+ 4^1  - k ( x 2 + y2)j2#2^>. 


( 16) 


This  energy  density  is  invariant  under  rotations  about  the 
z— axis,  and  the  Roynting  vector  has  a component  in  the  positive 
z— direction  and  one  radial  component  (perpendicular  to  the 
z— axis).  For  a light  pulse,  $ has  the  form 


= t (z)sin (kz) , 


$ (z  ) 


(17) 


where  k is  the  wave  number  o-f  the  light  produced  by  a laser. 

3 . Computer  Programs 

A computer  program  named  BALL  allows  us  to  calculate  the 
fields  by  performing  the  integrations  in  eqs.  (9)  and  (10).  This 
program  calls  a subroutine  ICS  where  the  unit  sphere  is  covered 
by  the  appropriate  number  of  patches  required  to  carry  out  the 
integrations,  and  a subroutine  POINT  that  is  used  to  step  through 
the  field  points.  The  output  of  BALL  is  stored  in  a file,  that 
is  then  used  as  input  to  separate  programs  that  present  the 
results  in  graphic  form.  These  programs  are  BALPLT,  which 
produces  a plot  of  the  energy  density  as  a function  of  the  field 
points,  and  8ALDSK,  which  prepares  a binary  interpolated  file 
that  can  be  shown  on  the  monitor  of  our  group's  image  processing 
facility  or  photographed  by  a special  computer-dr i ven  camera.  An 
additional  program,  BALCHK,  can  be  used  to  perform  a check  on  the 
accuracy  of  the  numerical  calculations  by  means  of  the 
conservation  of  energy  and  momentum.  These  programs  are  listed  in 
the  appendix. 

3. 1 Progr am  BALL 

We  use  this  program  to  do  the  integrations  required  to  find 

-4 

the  fields  at  a time  t and  at  a point  x for  a given  set  of 
initial  fields.  In  the  program  shown  in  the  appendix  the  pulse 
described  in  eqs.  (13)  and  (14)  has  a shape  given  by 

<|>(z)  = [exp  (az)  — exp  (£z  ) Isin  (kz  ) ® (—  z ) , (18) 
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i -f  a positive  k is  specified,  or 


$ ( z ) = zJ~Cexp(az)  — exp  (£z  ) 38  (—  z ) (19) 

if  k vanishes.  The  unit  step  function  6 makes  the  pulse  vanish 
for  positive  z at  t = 0.  For  the  pulsed  lasers  used  at  present, 
the  underlying  fields  oscillate  maybe  one  hundred  times  for  the 
duration  of  the  pulse,  but  this  number  is  being  reduced  in 
ongoing  experimental  research  programs.  Since  Maxwell's 
equations  in  free  space  are  invariant  under  changes  of  scale  in 
space  and  time,  the  same  calculations  can  be  applied  to  other 
el ectromagnet i c fields  that  may  be  true  pulses,  without  an 
underlying  sinusoidal  field. 

The  parameters  of  the  pulse,  the  region  over  which  the  field 
points  range,  and  other  parameters  of  the  calculation  are  given 
in  an  input  file  BALIN.XXX,  where  XXX  is  the  qualifier  that 

■4 

identifies  a particular  run.  The  initial  values  of  t and  x as 
well  as  the  increment  of  any  one  or  two  of  these  variables  are 
specified  in  this  file.  The  output  file  BALOUT. XXX  contains  the 
values  of  the  components  of  E and  cB. 

The  integration  is  performed  by  summing  the  contributions  to 
the  field  components  of  the  integrands  calculated  at  the  center 
of  each  patch  multiplied  by  the  area  of  the  patch. 

Different  versions  of  the  program  have  been  tried.  Some  o-f 
the  variations  change  the  way  the  patches  on  the  sphere  are 
determined  (a  simpler  subroutine  can  be  used  at  the  expense  o-f 
computing  time  or  accuracy) , start  from  other  initial  fields,  use 
double  precision  variables  to  do  the  computations  (improves  t vu* 
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accuracy  of  the  results  for  times  significantly  larger  than  the 
duration  of  the  initial  pulse),  or  save  only  the  energy  density 
instead  of  the  six  field  components  (saves  disk  space).  The 

program  contains  a provision  that  allows  us  to  restart  the 
calculation  without  much  loss  of  time  if  there  is  an 
i nterrupti on . 

3.2  Subroutine  F'O I.NJ 

This  subroutine  is  used  to  read  the  information  on  the 
region  of  space  and  time  where  the  fields  are  to  be  calculated, 

-4 

ana  to  generate  successive  values  of  x and  t. 

In  its  present  form,  the  program  requires  either  that  one 
variable  change  over  a given  range,  where  the  given  number  of 
points  is  to  be  distributed  uniformly,  or  that  two  variables 
change,  with  points  distributed  over  a rectangular  grid  uniformly 
in  each  direction.  A simple  variation  allows  the  program  to 
compute  the  fields  at  points  located  diagonally  in  the  xy— plane 
by  varying  both  x and  y simultaneously. 

3.3  Subroutine  IQS 

The  values  of  the  initial  fields  that  contribute  to  the 

■4 

fields  at  a point  x at  time  t are  all  on  the  Information— 

-4 

Collecting  Sphere  (ICS)  centered  at  the  field  point  x and  with  a 
radius  ct. 

If  the  initial  fields  have  the  form  given  by  eqs.  (13)  and 
(14),  with  a <£(z)  specified  in  (IS)  or  (19),  the  integrands  that 
contribute  significantly  to  the  result  are  confined  to  a finite 
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cylinder  whose  axis  is  the  z — axis,  whose  radius  is  a -function  of 

the  decay  constant  k of  the  gaussian,  and  which  is  limited  by  the 

planes  z = 0 and  z = z , where  z is  a function  of  the  decay 

m m 

constant  a of  the  double  exponential.  In  particular,  we  have  set 

2 

k ix^  + y*")  and  az  equal  to  a maximum  value  given  by  the  variable 
CUT  in  the  input  file.  We  thus  have  to  find  the  regions  of  the 
ICS  that  lie  inside  the  cylinder.  Depending  on  the  radius  of  the 
sphere  and  the  location  of  its  center,  the  whole  spherical 
surface  may  be  included,  one  or  two  regions  may  be  inside  the 
cylinder,  or  none  at  all.  Once  the  nature  and  the  limits  of 
these  regions  are  determined,  each  is  covered  by  patches 
determined  by  first  dividing  the  regions  into  strips  by  circles 
of  constant  polar  angle  Q and  then  dividing  each  strip  or  partial 
strip  into  patches  by  circles  of  constant  azimuthal  angle  0.  The 
sizes  of  these  patches  are  determined  in  two  different  ways.  A 
length  given  through  the  input  variable  DPATCH  is  used  to 
determine  the  magnitudes  of  the  sides  of  the  patches  except  when 
this  length  is  comparable  to  the  radius  of  the  sphere.  If  this 
is  the  case,  there  may  be  too  few  patches  to  assure  accurate 
values  for  the  integral,  and  a maximum  size  of  the  patches  is 
determined  by  dividing  the  sphere  into  N1  strips  and  dividing  the 
equatorial  strip  into  N2  patches,  N1  and  N2  being  values  supplied 
also  through  the  input  file. 

After  the  above  preliminary  computations  are  done  for  a 
given  time  and  field  point,  successive  calls  to  the  subroutine 
ICS  return  values  of  the  components  of  the  unit  vector  to  the 
center  of  a patch  and  the  solid  angle  subtended  by  the  patch. 
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given  by 


AO  = 2 a#  sin©  sin(?iA©).  (20) 

We  choose  a circular  cap  o-f  angular  diameter  equal  to  A©  as  the 
patch  at  either  of  the  poles?  the  correspond! ng  solid  angle  is 

Afl  = 4tt  sin^OiA©).  (21) 

n 

These  values  might  be  approximated  by  sin©A©Atf  and  n(?iA©)'w, 
respectively,  if  the  accuracy  needed  allows  this  substitution. 
We  nave  not  determined  precisely  when  this  is  allowable,  but  we 
found  that,  in  a particular  case,  the  switch  made  a significant 
difference. 

Subroutine  ICS  calls  the  subroutines  PATCH  to  determine  the 
components  of  the  unit  vector  for  a patch  on  a partial  strip, 
RPATCH  to  do  likewise  for  a full  strip,  and  FILL  to  transfer 
these  values  and  the  solid  angle  to  COMMON.  These  subroutines 
are  listed  together  with  ICS. 

3.4  Pr og^r am  BALPLT 


We 

use  this  program 

to  produce 

plots  of 

the  energy 

densi ty 

the 

el ectromagnetic 

field  as  a 

function 

of  the  one 

or  two 

variables  that  change  over  a given  interval.  The  input  comes 
from  the  files  BALIN.XXX  and  BALOUT.XXX,  and  the  output  is 
directed  to  a plotter  or  pr i nter /pi  otter . 

When  only  one  variable  changes,  the  data  on  the  dependent 
and  independent  variables,  as  well  as  labels,  are  passed  to  the 
subroutine  DRAW4,  a general  purpose  plotting  subprogram  described 
in  reference  C21.  When  two  variables  change,  the  number  of 
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points  in  either  variable  is  reduced  to  a maximum  of  100  if 


necessary,  and  then  the  data  are  passed  to  subroutine  PL0T3D, 
which  produces  three-dimensional  plots  and  is  described  in 
reference  [-33. 

The  nature  of  the  Green  function  for  the  wave  equation  is 
such  that  the  fields  propagate  exactly  with  the  speed  of  light. 
Thus,  the  field  in  the  xz  — plane  at  a time  t is  confined  to  a ring 
of  radius  ct  and  of  a width  equal  to  the  size  of  the  initial 
pulse.  When  the  grid  of  points  where  the  fields  are  calculated 
and  plotted  is  superimposed  on  this  ring,  the  appearance  of  the 
resulting  three-dimensional  plot  is  that  of  a large  number  of 
isolated  peaks  if  the  distance  between  the  points  is  comparable 
to  the  width  of  the  ring.  This  problem  can  be  minimized  by 
choosing  a density  of  points  that  is  high  enough. 

This  program  can  be  easily  adapted  to  show  all  components  of 
the  electromagnetic  field,  or  to  show  the  energy  density  when 
only  this  quantity  is  saved  in  the  file  BALOUT.XXX. 

3.5  Program  BALDSK 

This  program  was  developed  to  prepare  a raster  file  of 
values  proportional  to  the  energy  density  of  the  electromagnetic 
field  that  can  be  shown  on  the  monitor  of  the  image  processing 
facility  or  photographed.  We  assume  that  the  time  t is  fixed, 
and  that  the  grid  extends  over  a range  of  values  of  z and 
positive  values  of  x.  We  first  interpolate  in  the  x-direction  to 
256  points  and  then  reflect  this  set  of  values  about  the  z-axis 
to  take  advantage  of  the  symmetry  in  the  field,  and  then  w<> 


interpolate  to  512  such  records  in  the  z-direction.  The 
intensity  of  a single  graph  is  scaled  to  values  between  0 and 
255,  and  the  maximum  intensity  of  a particular  graph  in  a series 
for  varying  values  of  t is  scaled  to  fit  a range  between  150  and 
255 . 


3.6  Program  BALCHK 

To  use  BALCHK,  the  fields  have  to  be  computed  for  a fixed 
value  of  t and  for  values  of  z and  positive  x that  cover  the 
region  where  the  fields  are  significantly  different  from  zero. 
BALCHK  performs  a numerical  integration  to  compute  the  values  of 
the  total  energy  and  momentum  of  the  pulse;  these  values  allow  us 
to  check  the  validity  of  the  calculations  by  comparing  them  to 
the  initial  values  obtained  from  the  analytical  expressions  (15), 
(16),  and  (18)  or  (19).  The  x-  and  y-components  of  the  momentum 
vanish?  the  z-component  of  the  momentum  and  the  energy  density 
are  proportional  to 


dV  (E 


2 1 F 1 + 8(ec2  - 4 <x/3  + ff2)  + 1 1 

4k (a  + £)  5 40°’-'  ’ 


(22) 


J 


-4?  '?-+'? 

dV(E^  + c^B^) 

o o 


8 (g2  - 4ctff  + £2) 
(a  + 0)5 


64 

(a  + $ ) 5 


(23) 


12 


To  compute  these  integrals  numerically  at  later  times, 


we 


divide  the  space  where  we  compute  the  -fields  into  tori  o-f  square 
cross  sections  based  on  the  grid  in  the  xz  — plane,  and  we  multiply 
the  volume  o-f  each  torus  by  the  average  of  the  energy  density  or 
Poynting  vector  as  computed  at  the  four  corners  of  the  square  in 
the  grid. 

Since  all  the  contributions  to  the  energy  are  positive,  if 
the  grid  is  too  coarse  (that  is,  if  it  misses  much  of  the  ring 
where  the  fields  are  concentrated) , the  agreement  can  be  poor 
(that  is,  the  discrepancy  can  be  larger  than  a few  percent) 
without  the  values  of  the  field  being  incorrect.  In  these  cases, 
the  conservation  of  momentum  generally  shows  better  agreement 
than  the  conservation  of  energy. 

The  initial  values  of  the  z— component  of  the  Poynting  vector 
(15)  are  positive,  so  that  the  pulse  as  a whole  moves  in  the 
positive  z-direction. 

If  only  the  energy  density  of  the  field  is  saved,  we  can 
still  check  for  the  conservation  of  energy. 

4.  Sample  Outputs 

We  use  a pri nter /pi otter  to  obtain  the  outputs  from  the 
program  BALPLT.  In  fig.  1 we  show  three— di mensi onal  plots  of  the 

energy  density  of  the  simple  pulse  given  by  eq.  (19)  as 

-12 

functions  of  x and  z at  the  initial  time  and  after  1.2x10  s, 

-4 

which  corresponds  to  a distance  of  3.6x10  m.  The  pulse  moves 
mainly  forward  in  the  z— direction  and  to  the  sides,  and  the  peal- 

j 0 20 

intensity  decreases  from  3.3x10  to  5.9x10  in  arbitrary 
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uni ts. 


Figure  2 shows  the  same  features  for  a modulated  pulse 

specified  by  eq.  (18).  Although  the  wave  number  of  the 

modulation  is  not  much  bigger  than  the  size  of  the  pulse  (k  = 
5 —1  4 —1 

1x10  m and  a = 2x10  m ),  the  shape  of  this  pulse  changes 

little  over  the  chosen  period  of  time,  the  motion  is  mostly  in 

the  z-direction  and  the  maximum  intensity  decreases  from  0.74  to 

0.46  only.  The  superposition  of  the  platting  grid  on  the  ridges 

that  belong  to  the  energy  density  function  produces  additional 

peaks  in  the  plot.  Figure  3 shows  the  x— , y— , and  z-components 

of  the  fields  E and  cB  for  the  unmodulated  pulse.  They  are  shown 

-4  -4 

at  the  point  (1x10  , 1x10  , 0)  as  functions  of  time  between  0 

-12  -12 
and  1.2x10  s,  and  the  same  fields  at  the  time  1.2x10  s for 

-4 

y = 0 and  z = 2x10  as  a function  of  x between  the  values 
-4  -4 

—6x10  and  6x10  m (this  plot  exhibits  the  symmetry  about 

the  z— axis:  and  B^  are  antisymmetric  while  B_,  is  symmetric). 

Figure  4 shows  two  photographs  taken  with  a special  camera 
of  an  image  produced  from  a finer  grained  computation  of  the 
energy  density  shown  in  fig.  1,  with  the  fields  calculated  on  a 
256x256  grid.  This  camera  can  produce  a moving  picture  by 
photographi ng  a sequence  of  frames. 
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Fig,  1.  Energy  density  of  the  simple  pulse  as  a function  of 

^ 2 

and  z at  the  initial  time  and  after  1.2x10  s. 

cylindrical  symmetry  about  the  z-axis.  The  scale  of  distances 
shown  in  meters  and  the  maximum  intensities  in  arbitrary  units. 


x 

1 s 
1 s 
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■fig.  1 for  a modulated 


1 6 


Fig. 


Same 


as 


pul se. 


CO 


<s> 

ZT 


Fig.  3.  Plots  of  the  components  of  the  electromagnetic 
fields?  the  curves  are  labeled  from  1 to  6 to  identify  E^,  E.,, 

E_.,  and  B-.,  respectively.  Top!  functions  of  time  at  the 

point  (1x10  1x10  0)  . Bottom:  functions  of  x at  a time 

— 12  —4 

t = 1.2x10  for  fixed  y = 0 and  z = 2x10  m. 
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Fig  4.  Photographs  of  energy  densities  shown  in  fig.  1. 
The  pulse  is  propagating  mainly  to  the  left.  Symmetry  is  used  to 
extend  the  graph  to  negative  values  of  x. 


IQ 
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In  this  appendix  we  provide  listings  of  the  programs  and 
subroutines  BALL,  POINT,  ICS,  BALPLT,  BALDSK,  and  BALCHK. 

These  programs  use  some  general  purpose  software  provided 
with  the  Interdata  computer  system  or  developed  for  plotting.  We 
do  not  discuss  or  list  these  subroutines  here?  they  are  SYSIO, 
PLOT,  SYMBOL,  DRAW4,  and  PL0T3D. 
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PROGRAM 


BALL 


C 

C MICRO  AND  OPTICAL  METROLOGY  GROUP 

C NATIONAL  BUREAU  OF  STANDARDS 

C 

C EGON  MARX  11/29/83 

C 

C THIS  PROGRAM  IS  USED  TO  COMPUTE  THE  FIELDS  OF  A SPATIALLY  LOCALIZED 

C PULSE  AT  A TIME  T AND  AT  PRESCRIBED  FIELD  POINTS  WHEN  THE 

C INITIAL  FIELDS  ARE  GIVEN  AS  SOLENOIDAL  FUNCTIONS  OF  A VECTOR 

C VARIABLE  (X»  Y»  Z). 

C 

C R 1 » R2,  R3  = COMPONENTS  OF  THE  UNIT  VECTOR  TO  A PATCH 

C SOLANG  = SOLID  ANGLES  SUBTENDED  BY  A PATCH 

C ELI , EL2,  EL3  = ARRAYS  TO  STORE  COMPONENTS  OF  THE  ELECTRIC  FIELD 

C CBM1, CBM2, CBM3  = ARRAYS  TO  STORE  COMPONENTS  OF  THE  MAGNETIC  FIELD 

C CAPPA  = DECAY  CONSTANT  FOR  GAUSSIAN  IN  X-  AND  Y-DIRECTIONS 

C CUT  = CUTOFF  PARAMETER  FOR  INTEGRAND 

C ALPHA  = RISE  CONSTANT  FOR  DOUBLE  EXPONENTIAL  IN  Z-DIRECTION 

C BETA  = DECAY  CONSTANT  FOR  DOUBLE  EXPONENTIAL  IN  Z-DIRECTION 

C K = (REAL)  WAVE  NUMBER  FOR  MODULATION  IN  Z-DIRECTION 

C (NO  MODULATION  FOR  K = 0.  ) 

C DPATCH  = APPROXIMATE  SIZE  OF  PATCH 

C NI,  NJ  = NUMBER  OF  POINTS  IN  GRID  FOR  CALCULATION 

C 

DIMENSION  ELI (512), EL2< 512) , EL3( 512) , CBM1 ( 512) , CBM2<512>,  CBM3(512) 
LOGICAL  EX 
REAL  K, KZP 

CHARACTER  QAL*3, FL1*9, FL2#10 
COMMON  /RCOM/  R, R 1 , R2, R3, SOLANG 
COMMON  /ICOM/  N1,N2,  NX 

COMMON  /PARAM/  ALPHA, BETA, CAPPA, CUT, DPATCH, P I , TOP  I , P ITO, FOP  I , C 
INITIALIZE  VALUES 


NR=0 
C=3.  E8 

P IT0=2.  *ATAN< 1.  ) 

P 1=2. *PITO 
TOP  1=4.  #P ITO 
FOP  1=8.  *PITO 
FOP IM1=.  5/TOPI 

READ  IN  VALUES  OF  CONSTANTS  FOR  THE  PULSE 

OPEN ( 1 , FILE= 'BALL.  QAL 7 , STATUS= 'OLD  7 ) 

READ  (1,3)  QAL,  IFL 

FL1='BALIN.  7 //QAL 

FL2= 'BALOUT.  7 //QAL 

OPEN ( 3, F I LE=FL 1 , STATUS=  7 OLD  7 ) 

READ ( 3,  1 ) CAPPA, CUT, ALPHA,  BETA,  K,  DPATCH 
READ ( 3, 2 ) NI , NJ 

IF(NI.  GT.  512.  OR.  NJ.  GT.  512)  CALL  EXIT(l) 
NJ4=NJ*4 
N J44=NJ4+4 
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READ  FILE  QUALIFIER  AND  CHECK  FOR  RESTART 

I NQU I RE ( F I LE=FL2 , EXIST=EX) 

IF(  IFL.  EQ.  1.  AND.  EX)  THEN 

OPEN  < 2* FILE=FL2, ST ATUS=  7 OLD ' , ACCESS= 'DIRECT ' , 

1 FQRM= 'UNFORMATTED ', RECL=NJ4> 

INQUIRE ( 2, SI ZE=NR  > 

NX=NR/6 

IF(NX.  GE.  NI)  THEN 
CLOSE ( 2 ) 

STOP 

ENDIF 

ELSE 

OPEN  ( 2,  FILE=FL2,  STATUS=  'RENEW ' , RECL=NJ4,  BLOCKS  I ZE=N J44, 
1 FORM= 'UNFORMATTED',  ACCESS= 'DIRECT ' ) 

ENDIF 
CLOSE ( 2 ) 

GET  INITIAL  PULSE  DATA 

CALL  POINT < DUM,  DUM,  DUM,  DUM,  NI,  NJ,  . TRUE.  ) 

COMPUTE  THE  FIELDS  FOR  EACH  TIME  AT  THE  GIVEN  POINT 

DO  350  I=NX+1 , NI 
DO  300  J=l, NJ 

FIND  THE  COORDINATES  OF  THE  FIELD  POINT  AND  THE  TIME 

CALL  POINT<X,  Y,  Z,  T,  NI,  NJ,  . FALSE.  ) 

IF ( T.  NE.  0.  . OR.  Z.  GE.  0.  ) GO  TO  110 

COMPUTATIONS  FOR  T=0 


XP=X 

YP=Y 

ZP=Z 

XP2=XP**2 

YP2=YP**2 

ZP2=ZP**2 

CXP2=CAPPA*XP2 

CYP2=CAPPA*YP2 

CX2Y2=CXP2+CYP2 

AZP=ALPHA*ZP 

IF < AZP.  LT.  -25.  . OR.  CX2Y2.  GT.  25.  ) GO  TO  110 
EAZ=EXP ( AZP ) 

EBZ=EXP  < BETA-w-ZP ) 

DEXZ=EAZ— EBZ 

DEXZP=ALPHA*EAZ-BETA*EBZ 
IF < K.  EQ.  O.  ) THEN 
FZ=ZP2*DEXZ 

FPZ=ZP2*DEXZP+2.  *ZP*DEXZ 
ELSE 

KZP=K*ZP 


22, 


ooo  ooo  ooo  oooo»-ooo 


SKZ=SIN  <KZP ) 

CKZ=COS(KZP) 

FZ=DEXZ*SKZ 

FPZ=DEXZP*SKZ+DEXZ*K*CKZ 
END  IF 

TCXYP=2.  *CAPPA*XP*YP 
EC  X2Y2=EXP  < -C X2Y2 ) 

E1=YP*ECX2Y2*FPZ 

E2=-XP*ECX2Y2*FPZ 

E3=0. 

CB1=XP*ECX2Y2*FPZ 
CB2=YP*ECX2Y2*FPZ 
CB3=2. *<CX2Y2-1.  )*ECX2Y2*FZ 
GO  TO  210 
110  E1=0. 

E2=0 
E3=0. 

CB 1=0. 

CB2=0. 

CB3=0. 

IF ( T.  EQ.  0.  ) GO  TO  210 

ADD  THE  CONTRIBUTIONS  FROM  EACH  PATCH  ON  THE  UNIT  SPHERE 

CALL  ICS<X,  Y,  1,  T,  *200, *210) 

COMPUTE  THE  COMPONENTS  OF  THE  VECTOR 
FROM  THE  SOURCE  POINT  TO  THE  FIELD  POINT 

RU1=— R 1 
RU2=-R2 
RU3=-R3 
SLJ=SOLANG 
RV1=R*RU1 
RV2=R*RU2 
RV3=R*RU3 

COMPUTE  THE  COORDINATES  OF  THE  SOURCE  POINT 


XP=X— RV1 
YP=Y-RV2 
ZP=Z-RV3 

COMPUTE  THE  INTEGRANDS 

XP2=XP**2 

YP2=YP**2 

ZP2=ZP**2 

CXP2=CAPPA*XP2 

CYP2=CAPPA*YP2 

CX2Y2=CXP2+CYP2 

AZP=ALPHA*ZP 

FIND  DOUBLE  EXPONENTIAL  AND  ITS  DERIVATIVES 
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EAZ»EXP<AZP> 

EBZ«EXP(BETA*ZP> 

DEXZ«EAZ-EBZ 

DEXZP«ALPHA*EAZ~BETA*EB  z 
DEXZPP=ALPHA**2*EAZ-BETA**2*EBZ 
IF(K.  EQ.  0.  ) THEN 
FZ«ZP2*DEXZ 

FPZ«ZP2*DEXZP+2.  #ZP*DEXZ 
FPPZ»ZP2*DEXZPP+4.  *ZP*DEXZP+2.  *DEXZ 
ELSE 

KZP=K*ZP 

SKZ»SIN<KZP> 

CKZ»C0S(KZP) 

FZ=DEXZ*SKZ 

FPZ=DEXZP*SKZ+DEXZ*K*CKZ 

FPPZ=<DEXZPP-K**2*DEXZ)*SKZ+2.  *DEXZP*K*CKZ 
END  IF 

COMPUTE  FIELDS  AT  SOURCE  POINT 

TCXYP*2.  *CAPPA*XP*YP 
EC  X2Y2-EXP ( -CX2Y2 ) 

E01®YP*ECX2Y2*FPZ 

E02=-XP*ECX2Y2*FPZ 

E03«0. 

CB01=XP*ECX2Y2*FPZ 
CB02«YP*ECX2Y2*FPZ 
CB03=2.  * ( C X2Y2-1 . )*ECX2Y2*FZ 

COMPUTE  GRADIENTS  OF  FIELDS  AT  SOURCE  POINT 

DE01 1=-ECX2Y2*TCXYP*FPZ 
DE02 1 =EC X2Y2* ( 1 . -2.  *CYP2)*FPZ 
DE03 1 »EC  X2Y2*YP*FPP Z 
DE012=-ECX2Y2*( 1.  -2.  *CXP2)*FPZ 
DE022=-DE01 1 
DE032=-ECX2Y2*XP*FPPZ 
DE013=0. 

DE023-0. 

DE033®0. 

DB011«-DE012 
DB02 1 “-DE022 
DB03 1 s=-DE032 
DB012~DB021 
DB022«DE021 
DB032«DE031 

DB013=aECX2Y2*2.  *CAPPA»XP*<2.  -CX2Y2>*FZ 
DB023«ECX2Y2*2.  *CAPPA*YP*<2.  -CX2Y2)*FZ 
DB033=sECX2Y2*2.  *<CX2Y2-1.  >*FPZ 

COMPUTE  DIRECTIONAL  DERIVATIVES  ALONG  R 

R DEO 1 =R  V 1 *DEO 1 1 +R V2*DE02 1 +R V3*DE03 1 
RDEQ2=RV1*DE012+RV2*DE022+RV3*DE032 
RDE03=RV1*DE013+RV2*DE023+RV3*DE033 
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R DB  0 1 =R  V 1 *DBO 1 1 +R V2*DB02 1 +R V3*DB03 1 
RDB02=R V 1 *DB012+RV2*DB022+RV3*DB032 
RDB03=R V 1 *DBO 1 3+R V2*DB023+R V3*DB033 


COMPUTE  VECTOR  PRODUCTS 

EVRO 1 =E02*RU3-EG3*RU2 
EVR02=E03*RU 1 -EG 1 *RU3 
EVR03=E0 1 *RU2-E02*RU 1 
CBVR01=CB02*RU3-CB03*RU2 
CBVR02=CB03*RU1-CB01*RU3 
CBVR03=CB01*RU2-CB02*RU1 
RDEVR 1=RDE02#RU3 
RDEVR2=-RDE0 1 *RU3 
RDEVR3=RDE01*RU2-RDE02*RU1 
RDBVR 1 =RDB02*RU3-RDB03*RU2 
RDBVR2=RDB03*RU1-RDB01*RU3 
RDBVR3=RDB01*RU2-RDBQ2*RU1 

ADD  CONTRIBUTIONS  TO  THE  FIELDS  AT  THE  FIELD  POINT 


00 


DEl=SLJ*<E01-RDE01+<2. *CBVR01- 

DE2=SL J# ( E02-RDE02+  < 2.  *CBVR02- 

DE3=SLJ# ( E03— RDE03+  <2.  *CBVR03- 

DB1=SLJ# ( CB01 — RDB01  — ( 2.  *EVR01- 

DB2=SL J# ( CB02-RDB02- ( 2.  *EVR02- 

DB3=SL J* ( CB03— RDB03— ( 2.  *EVR03- 

E1=E1+DE1 

E2-E2+DE2 

E3=E3+DE3 

CB1=CB1+DB1 

CB2=CB2+DB2 

CB3=CB3+DB3 

GO  TO  120 

E1=E1*F0PIM1 

E2=E2*F0PIM1 

E3=E3*F0PIM1 

CB1=CB1#F0PIM1 

CB2=CB2*F0PZM1 

CB3=CB3*F0PIM1 


■RDBVR  1 ) ) 
■RDBVR2)  ) 
■RDBVR3 ) ) 
■RDEVR1 ) ) 
■RDEVR2)  ) 
•RDEVR3)  ) 


SAVE  THE  COMPONENTS  OF  THE  FIELDS 

IF  MAGNITUDE  LESS  THAN  l.E-39,  SET  EQUAL  TO  ZERO  TO  AVOID 
UNDERFLOW  WHEN  SQUARING 

210  IF(ABS<E1 ).  LT.  1.  E-39)  E1=0. 

IF(ABS(E2).  LT.  1. E-39)  E2=0. 

IF(ABS<E3). LT.  1. E-39)  E3=0. 

IF<ABS(CB1 ).  LT.  1. E-39)  CB1=0. 

IF  ( ABS ( CBS ) . LT.  1.  E-39)  CB2=0. 

IF(ABS(CB3>.  LT.  1. E-39)  CB3=0. 

ELI < J)=E1 
EL2 ( J ) =E2 
EL3 ( J ) =E3 
CBM1 < J)=CB1 
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CBM2< J)=CB2 
CBM3< J)=CB3 


300  CONTINUE 

C 

C SAVE  ARRAYS 

C 


OPEN (2, FXLE=FL2, STATUS" 'OLD ' < ACCESS* "DIRECT 
1 FORM* ' UNFORMATTED 7 , RECL=NJ4 ) 

CALL  WRITUN(EL1, NJ,  NR+1 > 

CALL  WRITUN<EL2>  NJ<  NR+2) 

CALL  WRITUN(EL3,  NJ,  NR+3) 

CALL  WRITUN<CBM1,  NJ,  NR+4) 

CALL  WRITUN<CBM2,  NJ» NR+5) 

CALL  WRITUN<CBM3,  NJ,  NR+6) 

NR=NR+6 
CLOSE  < 2 ) 

350  CONTINUE 


1 

2 

3 


STOP 

FORMAT (E12.  5) 
FORMAT (414) 
FORMAT <C3,  /,  II  > 
END 


C 


SUBROUTINE  WR ITUN ( A,  N,  IR ) 


C 

C 

C 

C 


THIS  SUBROUTINE  IS  USED  TO  WRITE  UNFORMATTED  DATA  TO 
A DIRECT  ACCESS  FILE 


DIMENSION  A ( N > 
WRITE(2,  REC=IR)  A 
RETURN 
END 
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SUBROUTINE  POINTCX, Y, Z, T, NI, NU,  INF) 


C 

C MICRO  AND  OPTICAL  METROLOGY  GROUP 

C NATIONAL  BUREAU  OF  STANDARDS 

C 

C EGON  MARX  11/29/03 

C 

C THIS  SUBROUTINE  RETURNS  A SEQUENCE  OF  FIELD  POINTS 

C TO  THE  MAIN  PROGRAM  BALL. 

C 

LOGICAL  INP 

DATA  DX,  DY,  DZ,  DT , NT  Is  NT2/4*0.  , 2*0/ 

COMMON  /RCOM/  R, Rt, R2, R3, SOLANG 
COMMON  /ICON/  N1,N2,  NX 

COMMON  /PARAM/  ALPHA, BETA. CAPPA,  CUT,  BPATCH,  P I , TOP I , P ITO,  FOP  I , C 
GO  TO  200  TO  READ  NEW  INPUT  DATA  IF  "INP"  IS  TRUE. 

IF (INP)  GO  TO  200 
PASS  NEW  COORDINATES 


X = XX 
Y=YY 
Z=ZZ 
T=TT 

NJ1=NJ1+1 

C CHECK  FOR  END  OF  ROW 

IF ( NJ1 . LE.  NJ ) THEN 
C INCREMENT  VARIABLE 

I F ( NT 1 . EG.  1 ) THEN 
XX=XX+OX 

ELSE  IFCNT1.  EQ.  2)  THEN 

YY=YY+DY 

ELSE  IFCNT1.EQ.  3)  THEN 
ZZ— ZX+DZ 
ELSE 

TT“TT-8”DT 
END  IF 

ELSE 
IN  J 1 = 1 

C RESET  VARIABLE  TO  INITIAL  VALUE 

IF < NT 1 . EQ.  1 ) THEN 
XX=X  0 

ELSE  IF < NT 1.  EG.  2)  THEN 
YY=YQ 
ELSE 
ZZ=ZQ 
END  IF 

C INCREMENT  OTHER  VARIABLE 

IFCNT2.  EQ.  2)  THEN 
YY=YY+DY 

ELSE  XFCNT2.  EQ.  3)  THEN 
ZZ=ZZ+DZ 
ELSE 
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TT=TT+DT 
END  IF 
END  IF 
RETURN 

STATEMENTS  TO  READ  IN  FIELD-POINT  COORDINATES, 

THE  INITIAL  TIME,  NUMBER  OF  THETA  STRIPS, 

AND  MAXIMUM  NUMBER  OF  PATCHES  ON  STRIP 

200  READ<3, 1 > XX, YY, ZZ, TT 

C SAVE  INITIAL  VALUES 

XO=XX 
YO=YY 
Z0=ZZ 

READ ( 3, 3 ) Nl, N2 

C READ  PULSE  SIZE  AND  DURATION 

READ(3, 1 > DELX, DELY, DELZ, DELT 
NJ1  = 1 

C CHECK  FOR  PROPER  GRID  VALUES 

IF ( NI . LE.  0.  OR.  NJ.  LE.  1 ) CALL  EXIT(7) 

C COMPUTE  COORDINATE  INCREMENTS 

NN=NJ 

IF  (DELT.  NE.  0.  ) THEN 
DT=DELT / ( NN-1 ) 

NN=NI 
NT  1=4 
NT2=4 
END  IF 

IF  (DELZ.  NE.  0.  ) THEN 
DZ=DELZ/ (NN-1 ) 

NN=NI 
NT  1=3 

NT2=MAX0 ( NT2, 3) 

END  IF 

IF  (DELY.  NE.  0.  ) THEN 
DY=DELY/ (NN-1 ) 

NN=NI 
NT  1=2 

NT2=MAX0 ( NT2, 2) 

END  IF 

IF  (DELX.  NE.  0.  > THEN 
DX=DELX/(NN-1 ) 

NT1  = 1 

NT2=MAX0 ( NT2, 1 ) 

END  IF 

C STOP  IF  ALL  DELX  ARE  ZERO 

IF ( NT1 . EG.  0 ) THEN 

WRITE ( 6, * ) 'ALL  INCREMENTS  ARE  ZERO' 

CALL  EXIT ( 8 ) 

END  IF 

C GIVE  WARNING  IF  ONLY  ONE  INCREMENT  IS  GIVEN  BUT  NI  IS  NOT  1 

IF ( NT1 . EG.  NT2.  AND.  NI.  GT.  1 > WRITE<6,  *) 

1 'WARNING:  ONLY  ONE  INCREMENT  GREATER  THAN  ZERO' 

IF ( NX.  EG.  0 ) RETURN 
C 
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260 

1 

3 


ADVANCE  COORDINATES  IF  RESTARTING 

IF ( NI . EQ.  1 ) CALL  EXIT<9) 

DO  260  1=1, NX 
IF(NT2.  EQ.  2)  THEN 
YY=YY+DY 

ELSE  IF(NT2.  EQ.  3)  THEN 
ZZ=ZZ+DZ 
ELSE 

TT=TT+DT 
END  IF 
CONTINUE 
RETURN 

FORMAT ( E12.  5) 

FORMAT < 214) 

END 
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SUBROUTINE  ICS< X,  Y,  Z,  T,  *,  * ) 


C 

C MICRO  AND  OPTICAL  METROLOGY  GROUP 

C NATIONAL  BUREAU  OF  STANDARDS 

C 

C EGON  MARX  i 1/29/83 

C 

C THE  FOLLOWING  GROUP  OF  SUBROUTINES  IS  USED  TO  COMPUTE  A PATCH 

C COVERING  FOR  THE  INFORMATION-COLLECTING  SPHERE,  AND  THEY  RETURN 

C TO  THE  MAIN  PROGRAM  BALL  THE  VALUES  OF  THE  COMPONENTS  OF  THE 

C UNIT  VECTOR  TO  A PATCH  AND  THE  SOLID  ANGLE  SUBTENDED  BY  THE  PATCH 

C 

C DPATCH  = APPROXIMATE  LENGTH  OF  A SIDE  OF  A PATCH 

C N1  = NUMBER  OF  THETA  STRIPS  (ALTERNATIVE) 

C N2  = MAXIMUM  NUMBER  OF  ANGLES  PHI  (ALTERNATIVE) 

C IP  - POINTER  FOR  PATCHES 

C IRP  = POINTER  FOR  PATCHES 

C 

DIMENSION  PMIN ( 2 ) , PMAX ( 2 ) , PRP (2 ) , RPMIN(2) , RPMAX (2) 

COMMON  /RCOM/  R, Rl, R2, R3, SOLANG 
COMMON  /ICOM/  N1,N2,  NX 

COMMON  /PARAM/  ALPHA, BETA, CAPPA, CUT, DPATCH, PI,  TOPI,  PITO,  FOP I,  C 
DATA  LBL/O/ 

GO  TO  ( 150, 200, 250),  LBL 

IP=0 

IRP=0 

G=  1 80.  /PI 

C COMPUTE  RADIUS  OF  CYLINDER  WHERE  INITIAL  DATA  ARE  SIGNIFICANT 

RHO=SQRT ( CUT/CAPPA ) 

C COMPUTE  COORDINATE  OF  BOTTOM  OF  CYLINDER 

ZM*3— CUT /ALPHA 

C COMPUTE  RADIUS  OF  INFORMATION-COLLECTING  SPHERE 

R=C*T 

C CHECK  THAT  Z IS  SUCH  THAT  THERE  MAY  BE  AN  INTERSECTION 

IF  ( Z.  LT.  R.  AND.  Z.  GT.  ZM-R  ) GO  TO  70 
RETURN  2 

C COMPUTE  DISTANCE  OF  FIELD  POINT  TO  Z-AXIS 

70  RP=SGRT(X**2+Y**2) 

C COMPUTE  AZIMUTH  ANGLE  OF  FIELD  POINT 

IF  ( X.  EG.  0.  . AND.  Y.  EQ.  0.  ) THEN 
PHI0=0. 

ELSE 

PHI 0=ATAN2 ( -Y,  -X) 

END  IF 

C COMPUTE  ANGULAR  SIZE  OF  PATCHES 

DTH=BPATCH/R 
NP=P I /DTH 

C COMPARE  TO  MINIMUM  NUMBER  OF  STRIPS 

IF ( NP.  LT. N1 ) DTH=PI/N1 
C COMPUTE  PATCH  SIZE  IN  PHI-DIRECTION 

DL2=DPATCH 

IF  (2-H-NP.  LT.  N2)  DL2=T0P  I*R/N2 
DTH2=DTH*.  5 

C COMPUTE  THETA  FOT  TOP  AND  BOTTOM  OF  CYLINDER 

THETAT=0. 
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THETAB=P I 

IF ( ABS ( Z ) . LE.  R>  THETAT=ACOS < -Z/R ) +1 . E-5 
IF ( ABS < ZM— Z ) . LE.  R)  THETAB=ACOS < ( ZM-Z  > /R ) 

C COMPUTE  DISTANCE  FROM  FIELD  POINT  TO  CYLINDER 

RPR=ABS(RP-RHO) 

C CHECK  IF  FIELD  POINT  IS  OUTSIDE  CYLINDER 

IF(RP. GT. RHO)  THEN 
C CHECK  FOR  INTERSECTION 

IF(R  LE.  RPR)  THEN 
RETURN  2 
ELSE 

C COMPUTE  THETA  FOR  INTERSECTION 

THETA 1=ASI N < RPR/R ) 

THET  A2=P I -THET A 1 

C CHECK  FOR  SECOND  INTERSECTION 

IF(R. LE.  RP+RHO)  THEN 

C CHECK  POSITION  OF  TOP  AND  BOTTOM  OF  CYLINDER 

IF ( THETA2.  LE. THETAT.  OR.  THETA 1 . GE.  THETAB ) THEN 
RETURN  2 
ELSE 

C SAVE  RANGE  OF  THETA  FOR  PARTIAL  STRIP 

IP=IP+1 

PMIN( IP )=AMAX1 (THETA1,  THETAT) 

PMAX ( IP ) =AMIN1 < THETA2, THETAB) 

PRP  < IP )=RP 
END  IF 
ELSE 

C COMPUTE  ANGLES  FOR  SECOND  INTERSECTION 

THETA3=AS I N ( ( RP+RHO ) /R ) 

THET  A4=P I— THET A3 

C CHECK  POSITION  OF  TOP  AND  BOTTOM  OF  CYLINDER 

IF (THETA2.  LE.  THETAT.  OR.  THETA  1 . GE.  THETAB.  OR. 

1 (THETA3.  LE.  THETAT.  AND.  THETA4.  GE.  THETAB) ) THEN 

RETURN  2 
ELSE 

IF< THETAT.  LT.  THETA3)  THEN 
C SAVE  RANGE  OF  THETA  FOR  PARTIAL  STRIP 

IP=IP+1 

PMIN( IP >=AMAX1 (THETAT, THETA1 ) 

PMAX ( IP )=AMIN1 (THETAB,  THETA3 ) 

PRP ( IP  >=RP 

IF ( THETAB.  GT.  THETA4)  THEN 
C SAVE  RANGE  OF  THETA  FOR  PARTIAL  STRIP 

IP=IP+1 

PMIN(IP) =THETA4 

PMAX ( IP ) =AMIN1 (THETAB,  THETA2 ) 

PRP (IP )=RP 
END  IF 
ELSE 

C SAVE  RANGE  OF  THETA  FOR  PARTIAL  STRIP 

IP=IP+1 

PMIN ( IP >=AMAX1 (THETAT,  THETA4 ) 

PMAX ( IP) =AMIN1 (THETAB,  THETA2 ) 

PRP ( IP ) =RP 
END  IF 
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END  IF 
END  IF 
END  IF 

C IF  THE  FIELD  POINT  IS  INSIDE  THE  CYLINDER 

ELSE 

C CHECK  IF  CYLINDER  DOES  NOT  INTERSECT  SPHERE  IN  TWO  CURVES 

IF ( R.  LE.  RP+RHO)  THEN 

C CHECK  IF  SPHERE  IS  INSIDE  CYLINDER 

IF(R.LE.  RPR)  THEN 

C SAVE  RANGE  OF  THETA  FOR  FULL  STRIP 

IRP=IRP+1 
RPMIN<  IRP  )*=THETAT 
RPMAX ( IRP )=THETAB 

C IF  SPHERE  INTERSECTS  CYLINDER  IN  ONE  CURVE 

ELSE 

C FIND  EXTREME  ANGLES 

THETA 1=ASIN< RPR /R ) 

THET  A2-P I -THET A 1 

C CHECK  POSITION  OF  TOP  OF  CYLINDER 

IF< THETAT.  LT.  THETA 1 ) THEN 
C SAVE  RANGE  OF  THETA  FOR  FULL  STRIP 

IRP=IRP+1 
RPMIN< IRP  )=THETAT 
RPMAX ( IRP) =AMIN1 ( THETA 1 , THETAB ) 

END  IF 

IF  (THETAT.  LT.  THETA2.  AND.  THETAB.  GT.  THETA  1 > THEN 
C SAVE  RANGE  OF  THETA  FOR  PARTIAL  STRIP 

IP=IP+1 

PMIN< IP )=AMAX1 < THETA li  THETAT) 

PMAX( IP)=AMIN1 (THETA2,  THETAB) 

PRP ( IP ) — — RP 
END  IF 

IF ( THETAB.  GT.  THETA2)  THEN 
C SAVE  RANGE  OF  THETA  FOR  FULL  STRIP 

IRP=IRP+1 

RPMIN ( IRP)=AMAX1 (THETA2,  THETAT) 

RPMAX ( IRP) =THETAB 
END  IF 
END  IF 

C IF  SPHERE  INTERSECTS  CYLINDER  IN  TWO  CURVES 

ELSE 

C COMPUTE  TWO  EXTREME  ANGLES 

THETA3=ASIN< < RP+RHO ) /R ) 

THET  A4=P I —THET A3 

C CHECK  POSITION  OF  TOP  AND  BOTTOM  OF  CYLINDER 

IF (THETAT.  GE.  THETA3.  AND.  THETAB.  LE.  THETA4 ) THEN 
RETURN  2 
ELSE 

C COMPUTE  OTHER  TWO  ANGLES 

THETA 1=ASIN (RPR /R ) 

THET  A2=P I -THET A 1 
C CHECK  TOP  OF  CYLINDER 

IF (THETAT.  LT.  THETA 1 ) THEN 
C SAVE  RANGE  OF  THETA  FOR  FULL  STRIP 

IRP=IRP+1 
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RPMIN < IRP  ) s=THETAT 

RPMAX ( IRP )=AMIN1 (THETAB,  THETA1 ) 

END  IF 

IF ( THETAT.  LT.  THETA3.  AND.  THETAB.  GT. THETA1 ) THEN 
C SAVE  RANGE  OF  THETA  FOR  PARTIAL  STRIP 

IP=IP+1 

PMIN( IP )=AMAX1 (THETAT,  THETA 1 ) 

PMAX< IP )=AMIN1 (THETAB,  THETA3 ) 

PRP ( IP ) =-RP 
ENDIF 

IF (THETAT.  LT.  THETA2.  AND. THETAB.  GT.  THETA4 ) THEN 
C SAVE  RANGE  OF  THETA  FOR  PARTIAL  STRIP 

IP=IP+1 

PMIN( IP)=AMAX1 (THETAT,  THETA4) 

PMAX( IP >=AMIN1 (THETAB,  THETA2) 

PRP ( IP ) =-RP 
ENDIF 

IF(THETAB.  GT.  THETA2)  THEN 
C SAVE  RANGE  OF  THETA  FOR  FULL  STRIP 

IRP=IRP+1 

RPMIN( IRP >=AMAX1 (THETAT,  THETA2 ) 

RPMAX ( IRP )=THETAB 
ENDIF 
ENDIF 
ENDIF 
ENDIF 

C END  OF  GEOMETRICAL  CONSIDERATIONS  BETUEEN  SPHERE  AND  CYLINDER 

C CHECK  FOR  PARTIAL  STRIPS 

140  IF( IP. EQ.  0)  GO  TO  190 
C INDICATE  PROCESSING  OF  PARTIAL  STRIP 

LBL=1 

C GET  DATA  FOR  STRIP 

THMIN=PMIN( IP) 

THMAX=PMAX( IP) 

PR=PRP( IP) 

C DECREMENT  POINTER 

IP=IP-1 

C COMPUTE  COMPONENTS  OF  UNIT  VECTOR  TO  PATCH 

150  CALL  PATCH(THMIN,  THMAX,  DTH,  DL2,  RHO,  PR,  PHIO, *140) 

RETURN 

C CHECK  FOR  FULL  STRIPS 

190  IF ( IRP.  EQ.  0 ) GO  TO  210 
C INDICATE  PROCESSING  OF  FULL  STRIP 

LBL=2 

C GET  DATA  FOR  STRIP 

THMIN=RPMIN( IRP) 

THMAX=RPMAX( IRP) 

C DECREMENT  POINTER 

IRP=IRP— 1 

C COMPUTE  COMPONENTS  OF  UNIT  VECTOR  TO  PATCH 

200  CALL  RPATCH(THMIN, THMAX, DTH, DL2, *190) 

RETURN 

C INDICATE  NO  MORE  STRIPS 

210  LEL=3 

RETURN 
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C ALL  DONE  FOR  THIS  FIELD  POINT 

250  LBL=0 

RETURN  1 
END 
C 

SUBROUTINE  PATCH < THMIN, THMAX, DTH, DL2, RHQ,  RP,  PHIO,  *) 

THIS  SUBROUTINE  COMPUTES  THE  COMPONENTS  OF  THE  UNIT  VECTOR 
TO  A PATCH  IN  A PARTIAL  STRIP 

COMMON  /RCOM/  R, Rl, R2, R3, SOLANG 

COMMON  /PARAM/  ALPHA, BETA, CAPPA, CUT, DPATCH,  P I , TOP I , P ITQ,  FOP  I , C 
DATA  LBL/O/ 

GO  TO  (100,  110,  120), LBL 
C CHECK  FOR  FIELD  POINT  ON  AXIS  OF  CYLINDER  (NO  PARTIAL  STRIPS) 

IF ( RP.  EQ.  0)  RETURN  1 
C COMPUTE  RANGE  IN  THETA 

DELTH=THMAX-THMIN 
C COMPUTE  NUMBER  OF  STRIPS 

NN1=DELTH/DTH 
DT1=DELTH/ (NN1+1 ) 

DT2=DT1*.  5 
B = . 5/RP 

A = ( RH0**2-RP**2 ) *B 
TH=THM I N+DT2 

C COMPUTE  SOLID  ANGLE  SUBTENDED  BY  PATCH 

SLG=2. *DL2*SIN(DT2)/R 
100  ST=SIN(TH) 

RST=R*ST 

CT=COS(TH) 

C COMPUTE  SIZE  OF  PATCH  IN  AZIMUTHAL  DIRECTION 

0PHI=DL2/RST 

C COMPUTE  RANGE  IN  AZIMUTH  FOR  PARTIAL  STRIP 

DELPHI=ACQS ( B*RST-A/RST  > 

NN2=DELPHI /DPHI 
DPHI=2.  *DELPHI/(2*NN2+1 ) 

PHI 1=PHI0 
PHI2=PHI0 

C SAVE  COMPONENTS  OF  UNIT  VECTOR  AND  SIZE  OF  CENTER  PATCH 

CALL  FILL(ST, CT, PHIO,  SLG) 

1 = 1 

IF(NN2.  EQ.  0)  GO  TO  130 
C INDICATE  ADDITIONAL  PATCHES 

LBL=2 
RETURN 

C INCREMENT  ANGLES 

110  PHI 1=PHI 1+DPHI 

PHI2=PHI2— DPHI 

C SAVE  COMPONENTS  OF  UNIT  VECTOR  AND  SIZE  OF  PATCH 

CALL  FILL(ST, CT, PHI1,  SLG) 

C INDICATE  SECOND  ANGLE 

LBL=3 
RETURN 

C SAVE  COMPONENTS  OF  UNIT  VECTOR  AND  SIZE  OF  PATCH 

120  CALL  FILL(ST, CT, PHI2, SLG) 
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1 = 1 + 1 

C CHECK  IF  DONE  WITH  STRIP 

IFd.LE.  NN2)  THEN 

C INDICATE  BACK  TO  FIRST  ANGLE 

LBL=2 
RETURN 
ENDIF 

C INCREMENT  THETA 

130  TH=TH+DT 1 

C CHECK  IF  DONE 

IF(TH. LE.  THMAX)  THEN 
C INDICATE  START  OF  NEXT  STRIP 

LEL=1 
RETURN 
ENDIF 

C INDICATE  ALL  DONE 

LBL=0 
RETURN  1 
END 
C 

SUBROUTINE  RPATCH ( THMIN,  THMAX* DTH,  DL2,  *) 

THIS  SUBROUTINE  COMPUTES  THE  COMPONENTS  OF  THE  UNIT  VECTOR 
TO  A PATCH  IN  A FULL  STRIP 

COMMON  /RCOM/  R, R 1 , R2, R3, SOLANG 

COMMON  /PARAM/  ALPHA, BETA, CAPPA,  CUT,  DPATCH,  P I , TOP I , P ITO. FOP  I , C 
DATA  LEL/O/ 

GO  TO  (80, 90,  100,  110), LBL 
IF  ( THMIN.  EQ.  0.  ) THEN 
C PROCESS  CAP  AT  TOP  OF  SPHERE 

DT2=AMIN1(DTH*.  5, THMAX ) 

S1=F0PI*SIN(DT2*.  5)**2 
CALL  FILL(0.  , 1.  , 0.  , SI  ) 

IF ( DT2.  EQ. THMAX)  RETURN  1 
THMIN=DT2 
LBL=1 
RETURN 
ENDIF 

80  IF (THMAX.  EQ.  PI ) THEN 

C PROCESS  CAP  AT  BOTTOM  OF  SPHERE 

DT2=AMIN1 (DTH*.  5, PI-THMIN) 

Si=F0PI*SIN(DT2*.  5>**2 
CALL  FILL ( 0.  , -1.  , 0.  , SI ) 

THMAX=P I— DT2 

IF ( THMIN.  GE.  THMAX)  RETURN  1 
LBL=2 
RETURN 
ENDIF 

C COMPUTE  RANGE  IN  THETA 

90  DELTH=THMAX-THM I N 

NN1=DELTH/DTH 
DT1=DELTH/ (NN1+1 ) 

DT2=DT1*.  5 
TH=THMIN+DT2 
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ST=SIN<TH> 

RST=R*ST 
CT=COB  < TH ) 

C COMPUTE  NUMBER  OF  PATCHES  IN  STRIP 

NN2=MAX0< INT<T0PI*RST/DL2+1.  ),  3) 

C COMPUTE  SIZE  OF  PATCH  IN  AZIMUTHAL  DIRECTION 

DPHI=T0PI/NN2 
PHI—O. 

C COMPUTE  SOLID  ANGLE  SUBTENDED  BY  PATCH 

SLG=2.  *DPHI*ST*SIN(DT2) 

1 = 1 

C INDICATE  NEXT  PATCH 

LBL=4 

C SAVE  COMPONENTS  OF  UNIT  VECTOR  AND  SIZE  OF  PATCH 

110  CALL  FILL < ST, CT,  PHI , SLG ) 

C INCREMENT  ANGLE 

PHI=PHI+DPHI 
1 = 1 + 1 

IF ( I . LE.  NN2 ) RETURN 
TH=TH+DT1 

IF(TH.  LT.  THMAX ) THEN 

C INDICATE  END  OF  STRIP  PROCESSING 

LBL=3 
RETURN 
ENDIF 

C INDICATE  ALL  DONE 

LBL=0 
RETURN  1 
END 
C 

SUBROUTINE  FILL ( ST, CT, PHI , SLG ) 

THIS  SUBROUTINE  COMPUTES  AND  STORES  THE  COMPONENTS  OF  THE 
UNIT  VECTOR  TO  THE  PATCH,  AND  THE  SIZE  OF  THE  PATCH 

COMMON  /RCOM/  R, R 1 , R2, R3, SOLANG 
R1=ST*C0S<PHI > 

R2=ST*SIN<PHI ) 

R3=CT 

SOLANG=SLG 
RETURN 
END 


36 


o o o o oon  noo 


PROGRAM 
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C 

C MICRO  AND  OPTICAL  METROLOGY  GROUP 

C NATIONAL  BUREAU  OF  STANDARDS 

C 

C EGON  MARX  11/29/83 

C 

C THIS  PROGRAM  IS  USED  TO  MAKE  A THREE-DIMENSIONAL  PLOT 

OF  THE  ENERGY  DENSITY  FROM  THE  OUTPUT  OF  PROGRAM  BALL  IF  TWO 
COORDINATES  VARY,  OR  A PLOT  OF  EACH  COMPONENT  OF  THE  ELECTROMAGNETIC 
C FIELD  IF  ONLY  ONE  COORDINATE  VARIES. 

C 

DIMENSION  EB( 100, 100), XLAB<2, 4), VAR (512) 

DIMENSION  ELI (512),  EL2(512),  EL3(512),  CBM1 (512),  CBM2(512), CBM3(512) 
INTEGER  PRLINE(20) 

CHARACTER  QAL*3, FL1*10,  FL2*9 
EQU I VALENCE  ( EB , VAR ) 

DATA  XLAB/  7 X— COORD.  ',  'Y-COORD.  ',  'Z-COORD.  ',  ' TIME  '/ 

INITIALIZE  THE  PLOTTING  ROUTINES 

READ  ( 5,  ■«■ ) IANG1,  IANG2 
THETA=I ANG1 
PHI=I ANG2 
READ ( 5, 3)  QAL 
FL1  = 'BALOUT.  ' //QAL 
FL2= 'BALIN.  7 //QAL 
OPEN <4, FILE=FL2,  STATUS= 'OLD ' ) 

READ ( 4,  2 ) NI,NJ 

CHECK  FOR  VALID  VALUES  FOR  PLOTTING  ANGLES 

IF ( NI . GT.  1.  AND.  < IANG1.  EQ.  0.  OR.  IANG2.  EQ.  0)  ) THEN 

WRITE (6,*)  'SUPPLY  NONZERO  VALUES  FOR  THETA  AND  PHI' 

CALL  EXIT(l) 

END  IF 
NJ4=NJ*4 
ND I=(NI— i)/i 00+ 1 
NDJ=(NJ-1 )/ 100+1 
READ (4,1)  XO, YO, ZO,  TO 
READ (4, 4)  IDUM, IDUM 
READ (4,  1 ) DELX, DELY,  DELZ,  DELT 
OPEN (2, F I LE=FL 1 , ST ATUS= ' OLD ' ) 

SELECT  VALUES  FOR  PLOTTING 
READ  IN  DATA  FILES 

11=0 

DO  110  1=1, NI 

CALL  READUN (ELI,  NJ ) 

CALL  READUN ( EL2,  NJ) 

CALL  READUN ( EL3,  NJ) 

CALL  READUN ( CBM1,  NJ) 

CALL  READUN ( CBM2,  NJ) 

CALL  READUN ( CBM3,  NJ) 
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100 

110 


120 

C 

C 

C 


130 

140 

150 

160 

170 

180 

C 

C 

C 


IF ( NI . EQ.  1 ) GO  TO  200 

I F ( MOD  ( I - 1 , ND I > . NE.  0 ) GO  TO  110 

11=11+1 

JJ=0 

DO  100  J=l, NJ,  NDJ 
J J3*  J J+  i 

EB(II, JJ)=EL1 ( J ) **2+EL2 ( J ) **2+EL3 ( J ) **2+ 
1 CBM1 ( J) **2+CBM2  < J ) *#2+CBM3 ( J ) **2 

CONTINUE 
CONTINUE 
UMAX=0. 

DO  120  1 = 1,  II 
DO  120  J=l, JJ 

UMAX=AMAX 1 ( UMAX , EB ( I , J ) ) 

CONTINUE 


PREPARE  A THREE-DIMENSIONAL  PLOT 

CALL  PLOTS < 0, 0,  0) 

ENCODE (PRL I NE,  13)  UMAX 

CALL  SYMBOL ( . 5,  9.  2,  . 1,  PRLINE,  0.  , 36) 

IF(DELY.  NE.  0.  . AND.  DELT.  NE.  0.  > GO  TO  130 

IF ( DELZ.  NE.  0.  . AND.  DELT.  NE.  0.  ) GO  TO  140 

IF (DELY.  NE.  0.  . AND.  DELZ.  NE.  0.  ) GO  TO  150 

IF ( DELX.  NE.  0.  . AND.  DELZ.  NE.  0.  ) GO  TO  160 

IF ( DELX.  NE.  0.  . AND.  DELY.  NE.  0.  ) GO  TO  170 

CALL  SYMBOL  ( . 5,  10.  , . 15, 

1 'ENERGY  DENSITY  AS  A FUNCTION  OF  X AND  T',0.  ,39) 

ENCODE (PRLINE, 7)  XO, YO, ZO, TO, DELX, DELT 
GO  TO  180 

CALL  SYMBOL  (.  5,  10.  , . 15, 

1 'ENERGY  DENSITY  AS  A FUNCTION  OF  Y AND  T',0.  ,39) 
ENCODE (PRLINE, 8)  XO, YO, ZO, TO, DELY, DELT 
GO  TO  180 

CALL  SYMBOL  (.  5,  10.  , . 15, 

1 'ENERGY  DENSITY  AS  A FUNCTION  OF  Z AND  T',0.  ,39) 

ENCODE (PRLINE, 9)  XO, YO, ZO, TO, DELZ, DELT 
CALL  SYMBOL  ( . 5,  10.  , . 15, 

1 'ENERGY  DENSITY  AS  A FUNCTION  OF  Y AND  Z',0.  ,39) 

ENCODE (PRLINE, 10)  XO, YO, ZO, TO, DELY, DELZ 
GO  TO  180 

CALL  SYMBOL  ( . 5,  10.  , . 15, 

1 'ENERGY  DENSITY  AS  A FUNCTION  OF  X AND  Z',0.  ,39) 

ENCODE (PRLINE, 11)  XO, YO, ZO, TO, DELX, DELZ 
GO  TO  180 

CALL  SYMBOL  (.  5,  10.  , . 15, 

1 'ENERGY  DENSITY  AS  A FUNCTION  OF  X AND  Y',0.  ,39) 

ENCODE (PRLINE, 12)  XO, YO, ZO, TO, DELX, DELY 
CALL  SYMBOL ( . 5,  9.  5,  . 1,  PRLINE,  O.  , 75) 

CALL  PLQT3D(EB,  0.  , 0.  , 100,  100,  1,  1,  II,  1,  1,  JJ,  1,  THETA,  PHI ) 
CALL  PLOT ( 0.  , 0.  , 999) 

STOP 


IF  ONLY  ONE  COORDINATE  CHANGES,  PREPARE  MATERIAL  FOR  PLOT 
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200  IF ( DELX. NE. 0. ) THEN 

V=X0 

DV=DELX / ( N J- 1 ) 

NTV=1 
END  IF 

IF (DELV.  NE.  0.  ) THEN 
V=Y0 

DV=DELY/(NJ-1 ) 

NTY=2 
END  IF 

IF (DELZ.  NE.  0.  ) THEN 
V=Z0 

DV=DELZ/(NJ-1 ) 

NTY=3 
END  IF 

IF  ( DELT.  NE.  0.  ) THEN 
V=TO 

DV=DELT / < N J— 1 > 

NTY=4 
END  IF 

DO  210  1=1, NJ 
VAR ( I ) =V 
V=V+DV 
210  CONTINUE 

ENCODE (PRLINE, 6)  XO, YO, ZO, TO 

CALL  DRAW4( 1, 7, 2, 4,  16,  17, XLAB ( 1, NTY) , 'FIELD  AMPLITUDES  ' , 

1 'ELECTROMAGNETIC  FIELDS,  1 : El,  2: E2,  3:  E3,  4:CB1,  5: CB2, 

2 ', PRLINE) 


PREPARE  ONE  PLOT  WITH  ALL  OF  THE  6 FIELD  COMPONENTS  AS  A 
FUNCTION  OF  THE  INDEPENDENT  VARIABLE 


7 

8 

9 

10 


CALL  DRAW4 < 2,  7,  1,  NJ,  71  100,  VAR,  ELI,  0.  , 0.  ) 

CALL  DRAW4 ( 2,  7,  1,  NJ,  72  7,  110,  VAR,  EL2,  0.  , 0.  ) 

CALL  DRAW4<2,  7,  1,  NJ,  '3',  120,  VAR,  EL3,  0.  , 0.  ) 

CALL  DRAW4 (2,  7,  1,  NJ,  '4',  130,  VAR,  CBM1,  0.  , 0.  > 

CALL  DRAW4<2,  7,  1,  NJ,  7 5 7 , 140,  VAR,  CBM2,  0.  , 0.  ) 

CALL  DRAW4<2,  7,  1,  NJ,  76  7,  150,  VAR,  CBM3,  0.  , 0.  ) 

CALL  DRAW4 ( 3,  7,  0,  0,  2,  -NJ,  ELI,  EL2,  2.  , 2.  ) 

CALL  PLOT < 0.  , 0.  , 999) 

STOP 

FORMAT ( E12.  5) 

FORMAT (//////,  214) 

FORMAT (C3) 

FORMAT < 214) 

FORMAT ( 7 INITIAL  VALUES:  X0=,,1PE8.  1, 

1 7,  T0=',  E8.  1 ) 

FORMAT < 7X0=7,  1PES.  1,  7 Y0=7,E8.  1, 

1 7 DELX=  7»  ES.  1,  ' DELT=  ' , E8.  1 ) 

FORMAT  < 7XG=7,  1PE8.  1,  ' Y0= 7 , E8.  1, 

1 7 DELY=  7 , E8.  1 , 7 DELT=  7 , E8.  1 ) 

FORMAT < 7X0=7,  1PE8.  1,  7 Y0=7,EB.  1, 

1 7 DELZ=  7 * E8.  1 , 7 DELT=7,E8.  1) 

FORMAT ( 7XG=7,  1PE3.  1,  7 Y0=7,E8.  1, 

1 7 DELY=  7 , E8.  1 , 7 DELZ=  7 , E8.  1 ) 


Y0= 
Z0=7,  E8.  1, 
Z0=7,  E8.  1, 
Z0=7,  E8.  1, 
Z0=7,  E8.  1, 


E8.  1,  7,  Z0= 

T0=  7 , E8.  1. 
T0=7,  E8.  1, 
T0=  7 , E8.  1, 
T0=  7 , E8.  1, 


6:  CB3 


, E8.  1, 
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11  FORMAT  ( ' XO-  ' > 1 PE8.  1 , ' Y0-SE8.  It'  Z0-SE8.li'  TO- S ES.  1, 
1 ' DELX-SE8.lt'  DELZ-SE8.  1) 

12  FORMAT < 'XO-S  1PE8.  1,  ' Y0-SE8.li'  Z0-SE8.lt  ' T0-SE8.  1, 
1 ' DELX-SES.lt'  DELY-SE8.  1) 

13  FORMAT ( 'MAXIMUM  ENERGY  DENSITY  IS  ',  1PE10.  3) 

END 

C 

SUBROUTINE  READUN(A,  N> 

THIS  SUBROUTINE  IS  USED  TO  READ  UNFORMATTED  DATA 

DIMENSION  A(N) 

READ ( 2 ) A 
RETURN 
END 
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C MICRO  AND  OPTICAL  METROLOGY  GROUP 

C NATIONAL  BUREAU  OF  STANDARDS 

C 

C EGON  MARX  1 1/29/83 

C 

C THIS  PROGRAM  PREPARES  FILES  TO  BE  PROCESSED  BY  THE  I2S 

C FACILITY  TO  BE  PUT  ON  A MONITOR  OR  PHOTOGRAPHED. 

C THE  INPUT  IS  THE  FILE  PRODUCED  BY  THE  PROGRAM  BALL. 

THE  OUTPUT  IS  A 512X512  SET  OF  INTEGERS  THAT  REPRESENT  THE 

C SCALED  INTENSITY  OF  THE  PICTURES. 

C WHEN  A SERIES  OF  OUTPUT  FILES  IS  PROCESSED  SIMULTANEOUSLY, 

C THE  MAXIMUM  INTENSITY  OF  EACH  OUTPUT  IS  SCALED 

C 

DIMENSION  EB ( 512) , IEB ( 128) , EB1 ( 512) , EB2< 512) 

DIMENSION  ELI (512), EL2(512),  EL3( 512) , CBM1 ( 512 ) , CBM2< 512 ) , CBM3 ( 512 > 
DIMENSION  AMX(IOOO) 

LOGICAL  EX 
CHAR AC TER *12  FL 
DATA  IFL/O/ 

OPEN  (5,  FTLE=  'BALDSK.  INP  7,  STATUS31  'OLD  7 ) 

READ (5,*)  N!,N2 

IF(N1.  GT.  N2.  OR.  N2.  GE.  10Q0)  CALL  EXIT<254) 

IF ( N1 . EQ.  N2 ) GO  TO  265 
C GET  OLD  FILE  OF  MAXIMA,  IF  IT  EXISTS 

INQUIRE (FILE- "BALDSK.  SAV7,  EXIST=EX  > 

IF (EX)  THEN 

OPEN (3, FILE= 'BALDSK.  SAV7,  STATUS* 'OLD  7 > 

GO  TO  240 
END  IF 

C OR  CREATE  A NEW  ONE  IF  NOT 

OPEN (3, FILE* 'BALDSK.  SAV7,  3TATUS= 'NEW ', 

1 FORM* ' UNFORMATTED  7 , RECL=4,  BL0CKSIZE=8> 

C OPEN  A FILE  FOR  OUTPUT  OF  MAXIMA 

OPEN ( 1 , F I LE=  ' BALDSK.  MAX',  STATUS*  7 RENEW  7 ) 

AMMAX=G. 

ALMAX=1 . E20 

C GET  THE  MAXIMA  FROM  THE  DIFFERENT  FILES 

DO  200  1 = 1,  1000 

FL=  BALOUT,  7//IT0C (I-i,  IDUM) 

OPEN (2, FILE=FL, STATUS= 7OLD 7 , ERR=220, SHARE=7SR07) 

C FIND  THE  SIZE  OF  THE  FILES 

IFdFL.  EQ.  0)  THEN 
IFL=1 

INQUIRE (2, SIZE=NI,  RECL=NJ4) 

NI=NI/6 
NJ=NJ4/4-l 
END  IF 
AMAX=Q. 

DO  110  11=1, N I 

CALL  READUN(EL1,  NJ) 

CALL  READUN(EL2,  NJ) 

CALL  READUN(EL3,  NJ) 

CALL  READUN(CEM1,  NJ) 
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CALL  READUNCCBM2,  NJ) 

CALL  READUN(CBM3»  NJ) 

DO  100  J=1,NJ 

AMAX=AMAX 1 ( AMAX , EL 1 ( J > **2+EL2  < J ) **2+EL3 ( J ) **2+ 
1 CBtll  ( J)#*24-CBM2(  J)**2+CBM3<  J)**2) 

100  CONTINUE 

110  CONTINUE 

CL0SE<2) 

AMX ( I ) -AMAX 

IF (AMAX.  LE.  AMMAX)  GO  TO  120 

AMMAX=AMAX 

NN=I 

120  IF ( AMAX.  G£.  ALMAX)  GO  TO  200 

AL!iAX=AMAX 
MM— I 

200  CONTINUE 
1=1000 

220  WRITE  < 3 ) 1-1 

DO  230  J=l, 1-1 
WRITE ( 3 ) AMX(J) 

WRITE (1,4)  J— 1 , AMX ( J ) 

230  CONTINUE 

WRITE (1,5)  NN-1,  AMMAX 
WRITE (1*6)  MM-1,  ALMAX 
WRITE(3)  AMMAX 
WRITE (3)  ALMAX 
GO  TO  260 

C RETRIEVE  INFORMATION  FROM  SAVED  FILES 

240  READ ( 3 ) I 

DO  250  J=l,  I 
READ(3)  AMX(J) 

250  CONTINUE 

READ(3)  AMMAX 
READ (3)  ALMAX 

260  C0EF=105.  / ( AMMAX-ALMAX ) 

C CREATE  RASTER  FILES  FOR  PICTURES 

265  DO  700  NN=N1  + 1 # N2+1 

FL= 'BALQUT.  '//ITOC (NN-1, IDUM) 

OPEN  ( 2,  FILE=FL*  ST  ATUS=  ' OLD  ' , SHARE='SRO' ) 

C FIND  INFORMATION  ON  FILE  SIZE 

IF ( IFL.  EQ.  0 ) THEN 
IFL=1 

INQUIRE (2, SI ZE=NI » RECL=NJ4> 

NI=NI/6 

NJ=NJ4/4-l 

ENDIF 

FL=  ‘'BALDSK.  '//IT0C(NN-1, IDUM) 

OPEN (4,  FILE=FL, FORM= ' B I NARY ' » ST ATUS=  7 RENEW ' , 

1 BLOCKS I ZE=512 ) 

WRITE (6, 11)  NN-1 

C FIND  MAXIMUM  FOR  SINGLE  PICTURE 

IF ( N1 . EQ. N2 ) THEN 
AMAX=Q. 

DO  2S0  11=1, NI 

CALL  READUN ( EL 1 , N J ) 
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270 

280 

C 


C 

C 

300 

C 


310 

320 

C 

330 

C 

c 

340 


CALL  READUN <EL2,  NJ) 

CALL  READUN <EL3,  NJ) 

CALL  READUN < CBM 1,  NJ) 

CALL  READUN(CBM2,  NJ) 

CALL  READUN <CBM3»  NJ) 

DO  270  J= 1 , N J 

AMAX=AMAX1 <AMAX,  ELI ( J)**2+EL2< J)**2+EL3< J)**2+ 

1 CBM 1 < J ) **2+CBM2  < J ) **2+CBM3  < J ) **2 ) 

CONTINUE 
CONTINUE 
FAC=255.  /AMAX 
REWIND  2 
ELSE 

FAC= ( COEF* ( AMX ( NN ) -ALMAX ) + 1 50.  ) /AMX < NN ) 

ENDIF 

READ  IN  FIRST  RECORD 
CALL  READUN  < EL  1 , N J ) 

CALL  READUN ( EL2, NJ) 

CALL  READUN < EL3,  NJ) 

CALL  READUN ( CBM1,  NJ) 

CALL  READUN ( CBM2i  NJ) 

CALL  READUN ( CBM3,  NJ) 

LIMIT  THE  NUMBER  OF  POINTS  TO  256  FOR  SYMMETRY 
NJ1=MIN0  < NJ, 256) 

SCALE  INTENSITIES 
DO  300  J=1,NJ1 

EB  ( J ) —FAC* < EL 1 ( J ) **2+EL2  < J ) **2+EL3 ( J ) **2+ 

1 CBM1 ( J)**2+CBM2< J)**2+CBM3< J)**2) 

CONTINUE 

INTERPOLATE  FIRST  RECORD 
I J=257 
EB02=EB ( 1 ) 

DO  320  XI=2, NJ1 
EBQ 1=EB02 
EB02=EB (II) 

K 1 = < 5 1 2- 1 J ) / < N J 1 + 1 — 1 1 ) 

FC— 1.  /KI 
DO  310  K=0, Kl-1 

EB2( I J)=FC*(EB01*<K1-K)+EB02*K) 

XJ=XJ+1 
CONTINUE 
CONTINUE 
EB2< IJ)=EB02 

FILL  IN  THE  REFLECTED  VALUES  ON  RECORD 
DO  330  11=1, 256 

EB2 (II) =EB2( 513-1 I ) 

CONTINUE 

PROCEED  WITH  OTHER  RECORDS,  INTERPOLATING  IN  OTHER  DIRECTION 
KJ=1 

DO  430  1=2, N I 

SHIFT  RECORD  BETWEEN  BUFFERS 
DO  340  K=l, 512 
EB1 ( K ) =EB2 ( K ) 

CONTINUE 

IJ=257 
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350 

C 


360 

370 

C 


380 

C 


C 


390 

400 

C 

C 

C 


420 

430 

C 


READ  ANOTHER  RECORD 

CALL  READUN (ELI*  N J ) 

CALL  READUN ( EL2,  NJ  > 

CALL  READUN ( ELS,  NJ) 

CALL  READUN ( CBM1,  NJ) 

CALL  READUN ( CBM2,  NJ) 

CALL  READUN (CBM3,  NJ) 

SCALE  NEW  RECORD 
DO  350  J=1 , NJ1 

EB  ( J ) =FAC* ( EL 1 ( J ) **2+EL2  < J ) **2+EL3  < J ) #*2+ 

1 CBM 1 ( J ) **2+CBM2 ( J ) **2+CBM3  < J > **2 ) 

CONTINUE 

INTERPOLATE  NEW  RECORD 
EB02=EB ( 1 ) 

DO  370  11=2, NJ1 
EB01=EB02 
EB02=EB (II) 

Kl  = ( 512— I J) / (NJ1+1— I I ) 

FC=1.  /K1 
DO  360  K=0, Kl-1 

EB2( IJ)=FC*(EB01*(K1-K)+EB02*K> 

IJ=IJ+1 
CONTINUE 
CONTINUE 
EB2( I J)=EB02 

FILL  IN  REFLECTION  OF  NEW  RECORD 
DO  380  11=1,256 

EB2( II )=EB2( 513-1 I ) 

CONTINUE 

INTERPOLATE  BETWEEN  RECORDS 
Ll  = < 512— KJ )/(NI  + l — I ) 

FC=1 . /LI 
DO  420  L=0, Ll-1 
IK=1 

DO  400  11=1, 128 
IVAL=0 

PACK  VALUES  IN  BYTES 
DO  390  K=1 , 4 

I VAL= I VAL*256+ 1 NT ( FC* ( EB 1 < I K ) * ( L 1 -L ) +EB2 ( X K > *L  > > 
IK=IK+1 
CONTINUE 
IEB (II ) =XVAL 
CONTINUE 

DO  SPECIAL  UNFORMATTED  I/O  TO  WRITE  TO  RASTER  FILE 

CALL  SYSIO( IPB, 57, 4, IEB, 512, 0) 

KJ=KJ+1 

CONTINUE 

CONTINUE 

DO  LAST  RECORD 
I K=  1 

DO  500  11=1, 128 
IVAL=0 

DO  450  K=l,  4 
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IVAL=IVAL*256+INT(EB2< IK) ) 

IK=IK+1 

450  CONTINUE 

IEB (II ) =X UAL 
500  CONTINUE 

DO  SPECIAL  UNFORMATTED  I/O  TO  WRITE  TO  RASTER  FILE 

CALL  SYSIO< IPB, 57,  4,  IEB, 512, 0) 

CLOSE (2) 

00  CONTINUE 
STOP 

FORMAT ( ' FILE  ',13,',  MAX= ' , 1PE12.  5) 

FORMAT  < //,  ' FILE  ',  13,  LARGEST  MAX=',  1PE12.  5) 
FORMAT < //,  ' FILE  ',13,',  SMALLEST  MAX= ', 1PE12. 5 ) 

1 FORMAT( IX,  13) 

END 

SUBROUTINE  READUNCA,  N) 

THIS  SUBROUTINE  IS  USED  TO  READ  UNFORMATTED  DATA 

DIMENSION  A(N) 

READ < 2 ) A 
RETURN 
END 
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PROGRAM 


BALCHK 


MICRO  AND  OPTICAL  METROLOGY  GROUP 
NATIONAL  BUREAU  OF  STANDARDS 

EGON  MARX  11/29/83 

THIS  PROGRAM  CHECKS  THE  CONSERVATION  OF  ENERGY  AND  MOMENTUM 
OF  THE  ELECTROMAGNETIC  FIELD  BY  COMPARING  THESE  QUANTITIES 
AT  TIME  T WITH  THEIR  INITIAL  VALUES. 

DIMENSION  EBH512),  EB2<  512),  EXCB  1(512  >,EXCB2(  512) 

DIMENSION  ELI  (512),  EL2(512>,  EL3(512>,  CBMK512),  CBM2012),  CBM3012) 
REAL  K,  K2 

CHARACTER  FL1*12, FL2*12,  QAL*3 
DATA  IFL/O/ 

P 1=4.  *ATAN(  1.  ) 

TOP  1=2.  *PI 

OPEN  FILE  WHERE  THE  QUALIFIERS  OF  FILES  TO  BE  PROCESSED  ARE  STORED 
OPEN ( 1 # FILE* 'BALCHK.  QAL',  STATUS* 'OLD ' ) 

OPEN  FILE  TO  SAVE  THE  OUTPUT 

OPEN (3, FILE* ' BALCHK.  OUT ' » STATUS* 'RENEW ' ) 

READ  (1,5,  END=400 ) QAL 

READ  PULSE  PARAMETERS 

FL 1 = ' BALIN.  7 //QAL 

OPEN (4, FILE=FL1,  STATUS* 'OLD ' ) 

READ (4, 1 ) CAPPA, DUM, ALPHA, BETA, K, DUM 

READ(4,  2)  NI,NJ 

READ(4, 1 ) XO, YO, ZO, TO 

READ(4, 2)  I DUM, I DUM 

READ (4, 1 ) DELX, DELY, DELZ, DELT 

DX=DELX/(NJ-1 ) 

DZ=DELZ/ ( NI-1 ) 

IF ( IFL.  EQ.  0 ) THEN 
IFL=1 


COMPUTE  INITIAL  VALUES  OF  THE  ENERGY  AND  MOMENTUM  OF  THE 
FIELDS  FROM  THE  ANALYTIC  EXPRESSIONS. 

C2=CAPPA*#2 

A2=ALPHA**2 

B2=BETA**2 

AB=ALPHA*BETA 

APB*ALPHA+BETA 

IF ( K.  EQ.  O.  ) THEN 

ERG=TQPI*(.  25/C2* ( . 25/ALPHA**3+8.  *(A2-4.  *AB 

1 +B2)/APB**5+. 25/BETA**3> 

2 +.  375/CAPPA*( 1.  /ALPHA#-* 5—64.  /APB**5+1.  /BETA**5>  ) 

EXB=P I#. 25/C2*(.  25/ALPHA**3+8.  *(A2-4.  *AB 

1 +B2)/APB**5+. 25/BETA**3) 

ELSE 

K2=K**2 

APB2=APB**2 

ERG=T0PI*K2*(.  25/C2*(. 25/ALPHA— 4.  * ( AB+K2 ) / ( APB* ( APB2+4.  *K2)) 
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250 


260 


270 


300 

310 


+ . 25/BETA ) +1 . /CAPPA*(.  25/ ( ALPHA# ( A2+K2 ) ) 

2 -4.  / ( APB*  < APB2+4.  *K2 ) > +.  25/ ( BETA* ( B2+K2  > ) > ) 

EXB=. 25*PX*K2/C2*( . 25/ALPHA-4.  *<AB+K2)/(APB*<APB2+4.  *K2> ) 

1 +.  25/BETA) 

END  IF 

WR ITE < 3#  4 ) ERG#  EXB 
END  IF 

FL2='BAL0UT.  V/QAL 

READ  IN  FIELD  VALUES  AND  COMPUTE  ENERGY  DENSITY  AND  POYNTING  VECTOR 
OPEN (2,  FILE-FL2# STATUS- "OLD ' > 

DO  THE  NUMERICAL  INTEGRATION  TO  FIND  THE  ENERGY  AND  MOMENTUM 


ERG=0. 

EXB=0. 

CALL  READUN ( EL 1 # N J ) 

CALL  READUN ( EL2# NJ) 

CALL  READUN ( EL3# NJ) 

CALL  READUN ( CBM1 # NJ) 

CALL  READUN ( CBM2# NJ) 

CALL  READUN < CBM3, NJ) 

DO  250  J=l, NJ 

EB2 ( J ) =EL 1 ( J ) **2+EL2 ( J ) **2+EL3 ( J ) **2+ 

1 CBM1 ( J ) **2+CBM2 ( J ) **2+CBM3 ( J ) **2 

EXCB2  < J ) =EL1 < J ) *CBM2 ( J ) —EL2  < J ) *CBM1 (J) 

CONTINUE 
DO  310  1=1, N I - 1 
DO  260  J=l, NJ 
EB1 < J)=EB2( J) 

EXCB1 (J)=EXCB2<J> 

CONTINUE 

CALL  READUN  < EL 1 # N J > 

CALL  READUN (EL 2, NJ) 

CALL  READUN ( EL3, NJ) 

CALL  READUN ( CBM1 , NJ ) 

CALL  READUN ( CBM2# NJ) 

CALL  READUN ( CBM3, NJ) 

DO  270  J=l,  NJ 

EB2  < J ) =EL 1 < J ) **2+EL2 ( J ) **2+EL3 ( J ) **2+ 

1 CBM1 < J)**2+CBM2< J)**2+CBM3< J)**2 

EXCB2 ( J ) =EL1 ( J ) *CBM2 ( J ) -EL2 ( J ) *CBM1 < J ) 

CONTINUE 

X1=X0 

DO  300  J~1 # NJ— I 
X2=X1+DX 

SURF=X2**2— X 1**2 

ERG=ERG+.  25*  < EB 1 ( J)+EB2( J)+EB1 ( J+l )+EB2( J+l ) )*SURF 
EXB-EXB+. 25*(EXCB1 ( J)+EXCB2( J)+EXCB1 ( J+l )+EXCB2( J+l ) )*SURF 
X 1=X2 
CONTINUE 
CONTINUE 
ERG=ERG*P I *DZ 
EXB=EXB*PI*DZ 
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UIRITE<3,  3)  FL1 * ERG»  EXB 
GO  TO  100 
400  STOP 
1 FORMAT (E12.  S) 

FORMAT ( 214 ) 

FORMAT(  IX,  C12,  2X,  3E15.  8) 

FORMAT  < //,  19X,  'ENERGY',  8X,  'MOMENTUM',/,  ' INITIAL  VALUES',  2E15.  8) 
FORMAT (A3) 

END 

SUBROUTINE  READUN( A, N) 

THIS  SUBROUTINE  IS  USED  TO  READ  UNFORMATTED  DATA 

DIMENSION  A(N) 

READ ( 2 ) A 
RETURN 
END 
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